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Abstract —This paper analyzes the Model Predictive Control 
technique for estimated Fluid Catalytic Cracking Unit (FCCU). 
It presents different estimated models of FCCU such as; impulse 
response model, estimated ARX model, state space model and 
output error model. These estimated models are then used as a 
model in Model Predictive Control (MPC) design. The main 
advantages of MPC is that it can handle hard input and output 
constraints and it can be used for Multi Input Multi Output 
processes (MIMO) without increasing the complexity in control 
design. MATLAB/Simulink is used to estimate the different 
models of FCCU and simulate the results of the controller. The 
simulation results shows that state space model estimated using 
N4SID logarithm provides better result for both identification 
and control. 

Index Terms —FCCU, MPC, N4SID, Output Error Model, 
System Identification. 

I. INTRODUCTION 

System identification is a process to build a mathematical 
model of dynamical systems from observed input and output 
data. These methods are widely used in industry for over a 
decade. It is being used in process control, aerospace, 
automotive, disk drives and embedded systems to create 
system models for any systems. The advantage of system 
identification is that it is quick and applicable to almost all 
systems [Billings, 2013]. 

On another hand Model Predictive Control (MPC) is an 
advanced optimized control method that has been widely used 
in process industries over the last two decades. It is a form of 
control in which the model of the process being controlled is 
required. The controller uses the model of the process and the 
output measurements to calculate the current control actions 
and predict the future behavior of the processes. The control 
action is calculated by minimizing the cost function at each 
sampling instant. 

This paper analyzes the estimated models of Fluid Catalytic 
Cracking Unit (FCCU) such as; impulse response model, 
estimated ARX model, state space model and output error 
model. These estimated models are then used as a model in 
Model Predictive Control (MPC) design. The response of the 
controller for these models are shown in this paper. 
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II. System Identification 

The procedure to identify the dynamical model for any 
systems are as follows. First we have to select an input signal 
and apply it to the plant to collect an output data. Input signal 
can be impulse signal, step signal, Pseudo Random Binary 
Noise (PRBS), Generalized Binary Noise (GBN), multiple 
sinusoids, etc. After getting the input and output data, the 
model structure of the plant is specified. There are three 
common types of models in system identification: white-box 
identification model, black-box identification model and 
grey-box identification model. White-box identification 
model structure is based on the first principles such as 
Newton’s law. In White-box identification; the model 
structure is completely known and the model parameters are 
estimated from the measured data. In grey-box identification 
the model structure is partially known from the first principles 
and the rest is developed from the measured data. In 
black-box identification, the model structure and its 
parameters are completely unknown and they are estimated 
using observed input and output data. After deciding the 
model structure, the system identification algorithms is used 
to estimate the mathematical models of dynamical systems. 
After the mathematical models of dynamical systems are 
developed, the result is validated. If the estimated model is 
not good enough then other estimation methods can be tried. 

Input signals play very important role in system 
identification because it is the only way to check the behavior 
of the process and collect the output data. In general we 
cannot introduce any random input to the process that is being 
estimated. We have to select an input signal that will carry 
enough knowledge about the system and affect smoothly all 
the operating frequencies. This can be achieved by using 
many different excitation signals, such as impulse signals, 
step signals, Pseudo Random Binary Sequences (PRBS), 
Generalized Binary Noises (GBN) and multiple sinusoids. 
The choice among these input signals depends on the type of 
identification technique used and the priori knowledge of the 
system under the test. 

In practice, for any controlled process the effect of the 
noises on the system should be minimized. The Input signal is 
the only freedom that the user have to determine the 
signal-to-noise ratio. Thus the amplitude of the input signal 
should be large enough in order to improve the signal-to-noise 
ratio but it cannot be too large that the output run away from 
the equilibrium point [Nelles, 2011]. 

Pseudo Random Binary Sequence (RBS) is a two stage 
deterministic signal with a periodic sequence of length N that 
switches between two levels, e.g. +a and -a !To generate these 
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sequences, there are two possibilities: first possibility is to use 
a quadratic residue code method suggested by Godfrey 
[Godfrey, 1993] and the second possibility is to use feedback 
shift register [Godfrey, 1993], [Eykhoff, 1974]. 

PRBS signals have been used for nonparametric model 
identification, such as: frequency response estimation and 
correlation analysis. In process control white noise signal is 
harmful to the actuator because it over emphasize the high 
frequency band at the cost of the low and the middle 
frequency band. In process control low pass character signals 
are preferred. Such signals can be obtained by increasing the 
clock period of the signal or by filtering the PRBS signals. 
The PRBS signal is preferred for system identification, 
because it excites all frequencies equally and imitate white 
noise in deterministic signal [Nelles, 2011]. In this paper 
PRBS signals are applied to FCCU to obtain an input-output 
data. 


A. Impulse Response Model 

Impulse response model refers to the models estimated 
using impulse response signal as an input. It is defined as 

6(t)=0 for all t^O (1) 

f 6(0 dt = 1 (2) 

From (1) and (2), it is clear that an ideal impulse function is a 
function that is zero everywhere but at the origin, where it is 
infinitely high. If a unit impulse function S( is applied as an 

input to a linear time invariant system (FTI) and the impulse 
response is denoted by gi, then a time delay input signal 
S(t - will give a time delay output signal represented by 
Eft-- If the input signal is represented as 

S(t - Tju(t)dt = u(tJ 

then the output impulse response will be 

g(t - x)u(x)dt = y(tj 

If gi is known, then for an input signal 
corresponding output signal can be computed. 




( 3 ) 

( 4 ) 
ui , the 


B. ARX Model 

ARX model relates the current output yi to a finite 
number of past outputs y(t— and inputs u(t-k). It is 
represented in linear difference equation as 

y(0 + a L y(t - 1) + - + a^yCt - nj 
= b t u(t - 1} +- — + b Db u(t - n b j + e(t) *■ ’ 

The advantage of ARX model is that it is the most efficient 
polynomial estimation methods available. It is considered to 
be the simplest way to represent the dynamic process driven 
by an input in the presence of error and disturbances. 

ARX model can be estimated using least square method. 
Rewriting (5) in regression form gives 

yCtej = 0(f) T 0 (6) 

where 

0(t) ' = [u(t - iXuft - 2), ..., u(t — n b 3] (7) 

e = [b L ,b 2 . b nt ]- (8) 

hence the model parameter can be estimated using the least 
square estimation given as 

§ = (9) 


C. State Space Model 


The state space representation of dynamical system given as 
x(k -f lj = AxfiO F BuQO F Ke(k3 (10) 

y(kj = Cx(k) -f Du(k) -f e(k) (11) 

where x (is the state vector, y( is the system output, u ( 
is the system input, e( is the stochastic error, A, B, C, D and 
K are the system matrices. 

Estimating the parameters in state space model is 
considered to be simple because the states or the order is the 
only parameter which needed to be estimated. 

Subspace identification is a method used to estimate the 
state space matrices A, B, C and D from an input and output 
data. It was proposed by van Overschee and de Moor. Further 
development to this method was done by Farimore in 1990 in 
which he proposed canonical variate analysis (CVA) 
[Farimore, 1990]. 

In 1994, van Overschee and de Moor proposed new 
numerical algorithm N4SID which identifies the mixed 
deterministic-stochastic systems. 

In N4SID, the observability and controllability of the 
system is not needed to be known in advance since the state 
space matrices are not calculated in their canonical form but it 
is calculated as the full state space matrices so there is no 
problem of identifiability [Overschee and Moor, 1994]. 

Fet consider a state space model of combined deterministic 
stochastic system given as 

Xfc + l = Ax k + Bu k + W k (12) 

y k = Cx k + Dll}, + 1 :^ (13) 

where A, B, C and D are the state space matrices, ^ is the 
state noise with covariance matrices E [v k v k ] = and is 


the output measurement noise with covariance matrices 
E [WfcwJ ] = and E [v k wj ] = 


If the system is observable then Kalman filter can be designed 
for the system to estimate the state variables according to 


% + i = Ax k + Bu k + K(y k - Cx k - Du k ) (14) 
where is the steady state Kalman gain that can be solved 
using Ricatti equation. 

P = APA t - K(C t PC 4- lOK 1- 4- Q (15) 

Denoting 

Sfc = y k - Cx k — Du k (16) 

Substituting (16) in (14) gives (10) and (11). 

The system described by (12) and (13) can be represented in 
the predictor form as 

Xk +1 = A K x k 4- B k z k (17) 

y k = Cx v + Du k + 6k (18) 

where 

2k = [u£,y£U K = A - KC, B k = [B — KD, K] (19) 


From (17) and (18), an extended state space can be 
formulated as 


Fftk) = F f x k -h G f z f _ t (kj -f BfUfCk) -h e f (kj 
where is the future horizon and 
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(24) 


By iterating (17) and (18), the following is obtained 


K k LpZp (kj-|- 


(25) 


where 




(26) 

(27) 

(28) 


Lp — . 3 k A k B k .... A]r 3 
ZpGO = l4.- 1^,-1 - 4- P ] 

By substituting (25) into (20) gives 
Yf 00 = I/Lp^ p 00 + 

+& Zf -1 (k) + D+Uftk) + &*- ( k) 

In subspace identification methods following assumptions 
are made: the eigenvalues of A — ] lies inside the unit 
circle, (A- Eis observable, (A. [B, h is controllable and 
is a stationary, zero mean white noise with covariance 

E[w k wJ] = 

Let consider an input vector x(kj and output vectory(kj, the 
linear regression can be expressed as 

y (k) = 0x00 -f v(k') 

which can be written in matrix form as 

[y(l)y(2) ...yGV)l 
= 0[x(l) x(2) .. . x(N) + \ 

Where V is a noise vector. 

By minimizing the following function 

J = || Y - 0Xll p 

Where |-1 is the F th norm, the least square solution is given 
as 

8 = YX'CXX 1 )" 

The model prediction is then given as 

Y = 0X = YX T CXX t ^- 1 X = YHx 
The least square residual is given as 


(29) 


(31) 


(32) 


(33) 


Y = Y - Y = Y(l - IIx) = YI1 - 

X 

Based form (10) and (11), an extended state space model can 
be formulated as 


Yf = r,X k T HfUj + GjEf 

where/ is the future horizon and 
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(40) 


(41) 


(42) 


The Kalman state is estimated from past input and output data 
based on equation 25 as 

= LpZpCX) + AgX k _ p 

where 

K k-p ~~ /k -p K k-p-§-i ■■■ K k-p-i-N-i] 

From equations 41 and 35 

Yf = r f L p Z p -f H f U f + G f E f 

= H|p Z p + H f U f + G f E f 

where = IfLp 

Under open-loop conditions, £> is uncorrelated to U f so 

1 T J (43) 

—E f Uj-» Oas N-* to ' 

Furthermore, £+- is uncorrelated to Z u from the Kalman filter 
theory. Therefore 

1 t (44) 

- E f Zj -» 0 as N -» go 

In N4SID we have to eliminate first eliminate Uf by 

post-multiplying n — on (42). 

u f 

1 111 ( 45 ) 

Y f n - = Hfp Zp n- +■ Hf Uf n - +■ G f E f n- 
Uf Uf Uf Uf 

T 

^ 3 Q) Then the noise term is removed by multiplying Zp from the 


result of equation (44). 


Y f n-Z P T = H* Zp n - Zp T + G f E f Zp T 
- H4- Z P n S Z P 


(46) 


and 


Hfp =Y f n-zJ v . p 




T-.-1 


(47) 


N4SID performs Singular Value Decomposition (SVD) on 
%Z p = F f L p Z p = USV 7 « (48) 

where S n contains the n largest singular value and 


(34) choosesZf = U r S 


L/2 


(49) 


D. Output Error Model 

The output-error model is defined, as follows 
<;(t) + %(t - 1) + ™ + - n f ) 

= biuCt - 1) + ™ + b Db u(t - n b ) 

where ^(t) is an undisturbed output, y(t) is the oputput at t, 
f Ll . ... , and b L . \ b are the unknown parameter that 

needed to be estimated. Re writing (49) in compact form gives 
y(t 0) = +- e ft J ( 5 °) 


where 


Tif 


= T f k q k = 1 + f t q 1 + - + f„ f q Df 


(51) 


k=c 
nt 


B(q) = > .b k q k = b L q 1 + - + b„. q nt 


k=i 


n t- 


(52) 


The output error model is estimated using regularization 
method. The parameters are estimated by minimizing the 
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mean square error given by a sum of systematic error (bias) 
and random error (variance). 

MSE = I Bias! 2 + VCE (53) 

III. Model Predictive Control 

Model Predictive Control (MPC) is considered to be an 
advanced optimized control method that has been widely used 
in process industries over the last two decades. The strategy 
that Model Predictive Controllers (MPC) uses to calculate the 

control actions characterized as ! At the k th sampling instant, 
the values of the control action, , for the next ‘M’ sampling 
instants. {u(kj,u(k T 1).. ..... u(k + M - fare calculated. 

They are calculated by minimizing the difference between the 
predicted outputs and the reference trajectories over the next 
T’ sampling instants while satisfying the constraints. In 
MPC, the control horizon ‘M’ and the predicted horizon T’ 
are the tuning parameters. After calculating the control moves 
for M sampling instants. The controller will implement the 
first control move uO At the next sampling instant, k+1, 
the control moves are recalculated for the next M sampling 
instants, [k+1 to k+M], and the first control 
move u(k + is implemented. These steps are repeated at 
each sampling instants. 
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predicted residuals and the control moves. A quadratic 

objective function can be written as 

P 


J(k) = ^"jlyCk + ilk) - r(k + OIIqq) 

i = i 
M—1 



+ ^ llAufk + ilk) II 
i = E 


where represents the model predicted output, is the 
reference setpoint, j is the change in manipulated input 
from one sample time to the next, is a weight for the 
changes in the manipulated input, is a weight for the 
changes in the predicted output, and the subscripts indicate 
the sample time . Pis a prediction horizon and M is a control 
horizon. 

The controller minimizes the objective function given in 
(54) to obtain the control moves. The objective function is 
minimized at each sampling instant using quadratic equation 
solvers. 

There are many different quadratic solvers available such 
as; interior point, active set, augmented Lagrangian, 
conjugate gradient, KIWIK algorithm, etc. The details of 
these algorithms are not considered in this paper. 

In general the quadratic problem can be formulated as 


+ -u T Hu) 


(55) 


Subject to Ax < b 

where is the optimal solution which gives the minimum, H 
is the positive definite Hessian matrix, A is a matrix of linear 
constraint coefficients, and b and fare the vectors. The H and 
A matrices are constants matrices which are calculated during 
the initialization of the controller and b and f vectors are 
calculated at the beginning of each sampling instant. 

The quadratic solver tries to find the minimum of the 
function given in (55) which satisfies the constraints. 

In MPC calculation, it is assumed that all the states of the 
process are measurable but that’s not true. In most of the 
application it is impossible to measure all the states and the 
estimation of the states are required. Observer, Kalman filter 
and Extended Kalman filter can be used for the state 
estimation. 


Fig. 1MPC Strategy 

In MPC strategy, the dynamical model of the plant being 
controlled is required. It is needed to calculate the future 
predicted output. The dynamical model of the plant must be 
accurate thus it must capture the dynamics of the processes 
completely. If the dynamical model of the process is not 
accurate then the result of the controller might not give the 
desired performance of the process outputs. There are many 
different types of models which can be used in MPC, such as 
impulse response, step response, transfer function models and 
sate space models. 

To obtain the control law, objective function is needed. 
Different objective function leads to different algorithms. 
There are several different choices for objective functions. 
The first one is a standard least-squares or quadratic objective 
function. The objective function is a sum of squares of the 


IV. Fluid Catalytic Cracking Unit 

The fluid catalytic cracking unit (FCCU) is complex 
interactive process in in petroleum refining industries. It 
takes chains of hydrocarbons and breaks them into smaller 
ones which allows refineries to utilize their crude oil 
resources more efficiently. It uses an extremely hot catalyst to 
crack the hydrocarbons into smaller ones. FCCU processes 
present challenging control problem because it has very 
complex kinetics of both cracking and coke burning reactions 
and it has strong interaction between the reactor and 
regenerator. A schematic of FCCU is shown in fig. 2. 

The FCCU contains of two main parts: the reactor and the 
regenerator. In reactor, the reaction between the mixture of 
hydrocarbon vapors and catalyst takes place. This reaction 
breaks the long molecules of hydrocarbons into smaller one 
which leaves from the top of the reactor. Steam is supplied to 
remove hydrocarbons from the catalyst. The cracking reaction 
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produces carbon materials and un-cracked organic materials 
known as coke that reduces the catalyst activity. The catalyst 
is taken into regenerator where it is regenerated by burning off 
the deposited coke with air. The regenerated catalyst is then 
taken to the reactor to repeat the cycle. The combustion of the 
coke in the regenerator produces a heat absorbed by the 
regenerated catalyst. This absorbed heat provides the energy 
for the vaporization and cracking reactions that take place in 
the catalyst riser. 



Fig. 2 Schematic of FCCU 


The FCCU considered here is developed from [Skogested, 
1993]. It has two manipulated variable and two controlled 
variables. The manipulated variables are: .which is the flow 
rate of regenerated catalyst and jwhich is the flow rate of air 
to regenerator. The controlled variable are: Twhich is the 
riser outlet temperature and 7 which is the regenerator 

temperature. 

V. Results and Discussions 

PRBS signals applied to the model of FCCU developed in 
Simulink. The developed model is based on the work of 
Skogested, 1993. This model has two manipulated variables 
and two controlled variables. The input-output data is then 
used to estimate the different models of FCCU. 

The estimated impulse model is shown in fig. 3 and fig. 4. 
From those figures, it is clear that the impulse model was 
unable to estimate the correct model. Impulse model 
produced a fit of 4.214% forT r , e , and -555.1% fit for T rs . 


Measured and simulated model output 



Fig. 3 Riser Outlet temperature, Impulse Model 


x iq 4 Measured and simulated model output 
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Fig. 4 Regenerator Temperature, Impulse Model 

ARX model is shown in fig. 5 and fig. 6. This model 
produced a fit of 86.8% for T ro , and 85.5% fit for 7^. 


Measured and simulated model output 



Fig. 5 Riser Outlet temperature, ARX Model 

Measured and simulated model output 



Fig. 6 Regenerator Temperature, ARX Model 


State Space Model is shown in fig. 7 and fig. 8.This model 
produced a fit of 87.83% for 7^, and 84.68% fit for . 


Measured and simulated model output 



Fig. 7 Riser Outlet temperature, State Space Model 
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Measured and simulated model output 



Fig. 8 Regenerator Temperature, State Space Model 


Output error model is shown in fig. 9 and fig. 10. This 
model produced a fit of 84.65% for T ro , and 80.8% fit for 



Measured and simulated model output 



Time 

Fig. 9 Riser Outlet temperature, Output Error Model 


Measured and simulated model output 



Fig. 10 Regenerator Temperature, Output Error Model 

Since the estimated impulse model was unable to produce a 
correct model, it was not used in MPC design. In MPC ARX 
model, state space Model and output error model are used. 
The results of MPC design is shown in fig. 11 and fig. 12 and 
fig. 13. The results show that the state space model estimated 
using N4SID produce better response than ARX model and 
ARX model produces better response than output error model 


Plant Outputs 



Rant Outputs 



1500 


Time (seconds) 

Fig. 12MPC Design State Space Model 

Rant Outputs 



-500 



Fig. 13MPC Design Output Error Model 
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VI. Conclusion 

In this paper, system identification techniques are used to 
estimate the dynamical models of FCCU. Impulse model, 
ARX model, state space model and output error model are 
estimated and compared. The result of the simulation shows 
that an Impulse model was unable to give a correct fit, 
however; ARX model, state space model and output error 
model produced a good fit. The estimated models then used in 
MPC design to control the riser outlet temperature and 
regenerator temperature of FCCU unit. MPC design is 
simulated in MATLAB. The result of the simulation shows 
that the estimated state space model using N4SID algorithm 
gives better result than ARX model and ARX model gives 
better result than output error model. 
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